function b = rhs_bb(V,T,d,dof_map,desc,Re,caseNum)
%
%   obtain the right hand side of the weak form of stream-function form NS 
%   equation. With the right hand side is given by function rhs(). 
%

    dim = max(max(dof_map));
    b = zeros(dim,1);

%   the quadreture points
    [qw,qp] = quad_rule(13);
    qx = qp(:,1); qy = qp(:,2);

    %   finite element data on master elements
    [phi,phi_x,phi_y,phi_xx,phi_xy,phi_yy] = fe_bb(d,qx,qy,desc);
    
    % geometry informations
    x13 = V(T(:,1),1) - V(T(:,3),1);
    x23 = V(T(:,2),1) - V(T(:,3),1);
    y13 = V(T(:,1),2) - V(T(:,3),2);
    y23 = V(T(:,2),2) - V(T(:,3),2);
    jac = (x13.*y23 - x23.*y13);
           
    nt = size(T,1);
    for k=1:nt
        % gradient value for the basis on each element
        phix = (y23(k)*phi_x - y13(k)*phi_y)'/jac(k);
        phiy = (x13(k)*phi_y - x23(k)*phi_x)'/jac(k);
        
        px = V(T(k,3),1) + qx*x13(k) + qy*x23(k);
        py = V(T(k,3),2) + qx*y13(k) + qy*y23(k);
        
        [f1,f2] = rhs(px,py,Re,caseNum);
        
        loc_b = phiy*(qw.*f1) - phix*(qw.*f2);
        
        % 0.5 is the area of the master element
        b(dof_map(:,k)) = b(dof_map(:,k)) + loc_b*jac(k)*0.5; 
        
    end
end